%sample code
%crack test after element has experienced crack
%element 1 formulation
function PatchTestT3Crack2
global H e nu

clc
H = -1; %-1
e = 2;
nu = 0.3;
ft = 0.91;

C = e/(1-nu*nu)*...
    [   1 nu 0;
        nu 1 0;
        0 0 (1-nu)/2    ];

%element 1
B1 = [   -1 0 1 0 0 0;
        0 0 0 -1 0 1;
        0 -1 -1 1 1 0   ];
w1 = 0.5;
dJ1 = 1;
l1 = 2/3;
A1 = 0.5;
n1 = [-1 0]';
m1 = [0 -1]';

Ge1 = -l1/A1*[n1(1) 0;0 n1(2);n1(2) n1(1)];

Gc1 = [  1 0;
        0 0;
        0 1    ]; % satisfy trivial relation (41)
%check trivial relation for Gc
xi0 = [1 1.1]';
dhat = [xi0' 0 0 0 0]';
disp(['trivial condition 1 = ' num2str(norm(B1*dhat+Gc1*xi0))]);

Ke1 = w1*dJ1*B1'*C*B1;
Kdxi1 = w1*dJ1*B1'*C*Gc1;
Kxid1 = w1*dJ1*Ge1'*C*B1;

J1 = [n1 m1];
Cg1 = [ H 0;0 0 ]*J1';
Kxixi1 = w1*dJ1*Ge1'*C*Gc1 + l1*J1'*Cg1*J1;

Kt1 = Ke1 - Kdxi1*inv(Kxixi1)*Kxid1;
%%%%end of element 1%%%%

%element 2
B2 = [   0 0 1 0 -1 0;
        0 -1 0 0 0 1;
        -1 0 0 1 1 -1   ];
w2 = 0.5;
dJ2 = 1;
l2 = 2/3;
A2 = 0.5;
n2 = [1 0]';
m2 = [0 1]';

Ge2 = -l2/A2*[n2(1) 0;0 n2(2);n2(2) n2(1)];

Gc2 = [  -1 0;
        0 0;
        0 -1    ]; %satisfy relation (41)
%check trivial relation for Gc
xi0 = [1 1.1]';
dhat = [0 0 xi0' 0 0]';
disp(['trivial condition 2 = ' num2str(norm(B2*dhat+Gc2*xi0))]);

Ke2 = w2*dJ2*B2'*C*B2;
Kdxi2 = w2*dJ2*B2'*C*Gc2;
Kxid2 = w2*dJ2*Ge2'*C*B2;

J2 = [n2 m2];
Cg2 = [ H 0;0 0 ]*J2';
Kxixi2 = w2*dJ2*Ge2'*C*Gc2 + l2*J2'*Cg2*J2;
% inv(Kxixi2)

Kt2 = Ke2 - Kdxi2*inv(Kxixi2)*Kxid2;
%%%%end of element 2%%%%

%element connectivity
con1 = [1 2 3 4 5 6];
con2 = [1 2 5 6 7 8];

%boundary condition
free_dofs = [6 8];
prescribed_dofs = [3 5];
fext = [0 0]';
ui = 0.5;
u_p = [ui ui]';

%initial condition
K = zeros(8);
K(con1,con1) = K(con1,con1) + Ke1;
K(con2,con2) = K(con2,con2) + Ke2;
ui0  = 0.45;
u_p0 = [ui0 ui0]';
Kuu = K(free_dofs,free_dofs);
Kup = K(free_dofs,prescribed_dofs);
u0 = Kuu\(fext-Kup*u_p0);
usystem0 = zeros(8,1);
usystem0(free_dofs) = u0;
usystem0(prescribed_dofs) = u_p0;
u = usystem0;

%system tangent stiffness matrix
Kt = zeros(8);
Kt(con1,con1) = Kt(con1,con1) + Kt1;
Kt(con2,con2) = Kt(con2,con2) + Kt2;

%crack initialization
xi1 = [0 0]';
xi2 = [0 0]';

%trigger displacement control
u(prescribed_dofs) = u_p;

Du = zeros(8,1);
tol = 1e-6;
i = 1;
cnt = 0;
while 1
    disp(['step ' num2str(i) '-----------------']);
    
    ue1 = u(con1);
    epsilon1 = B1*ue1 + Gc1*xi1;
    sigma1 = C*epsilon1;
    u_mu1_n = xi1(1);
    t1 = ft+H*u_mu1_n;
    Tg1 = [max(0,t1) 0]';
    renh1 = -Ge1'*sigma1*w1*dJ1 - J1'*Tg1*l1;
    if t1 < 0
        Kxixi1 = w1*dJ1*Ge1'*C*Gc1;
        Kt1 = Ke1 - Kdxi1*inv(Kxixi1)*Kxid1;
    end
    ri1 = B1'*sigma1*w1*dJ1+Kdxi1*(Kxixi1\renh1);
    
    ue2 = u(con2);
    epsilon2 = B2*ue2 + Gc2*xi2;
    sigma2 = C*epsilon2;
    u_mu2_n = xi2(1);
    t2 = ft+H*u_mu2_n;
    Tg2 = [max(0,t2) 0]';
    renh2 = -Ge2'*sigma2*w2*dJ2 - J2'*Tg2*l2;
    if t2 < 0
        Kxixi2 = w2*dJ2*Ge2'*C*Gc2;
        Kt2 = Ke2 - Kdxi2*inv(Kxixi2)*Kxid2;
    end
    ri2 = B2'*sigma2*w2*dJ2+Kdxi2*(Kxixi2\renh2);
    
    Kt = zeros(8);
    Kt(con1,con1) = Kt(con1,con1) + Kt1;
    Kt(con2,con2) = Kt(con2,con2) + Kt2;

    ri = zeros(8,1);
    ri(con1) = ri(con1) + ri1;
    ri(con2) = ri(con2) + ri2;
    
    du = Kt(free_dofs,free_dofs)\(fext-ri(free_dofs));
    u(free_dofs) = u(free_dofs) + du;
    
    Du(free_dofs) = du;
    
    due1 = Du(con1);
    dxi1 = Kxixi1\(renh1-Kxid1*due1);
    due2 = Du(con2);
    dxi2 = Kxixi2\(renh2-Kxid2*due2);
    
    xi1 = xi1+dxi1;
    xi2 = xi2+dxi2;
    
    conv_criteria = norm(du)/norm(u);
    
%     u_mu1
%     u_mu2
    renh1
    renh2
    Tg1
    Tg2
    xi1
    xi2
    sigma1
    sigma2
    epsilon1
    epsilon2
    i = i+1;
    
    cnt = cnt+1;
    if(conv_criteria < tol) 
        disp(['converged after ' ...
                num2str(cnt) ' loops']);
        break;
    end
end
u
% 
% %perform another iteration
% disp(['start new iteration']);
% i = 1;
% ui = 0.95;
% u_p = [ui ui]';
% %trigger displacement control
% u(prescribed_dofs) = u_p;
% cnt = 0;
% while 1
%     disp(['step ' num2str(i) '-----------------']);
%     
%     ue1 = u(con1);
%     epsilon1 = B1*ue1 + Gc1*xi1;
%     sigma1 = C*epsilon1;
%     xi1
%     u_mu1 = J1*xi1
%     Tg1 = [max(0,ft+H*u_mu1(1));0];
%     renh1 = -Ge1'*sigma1*w1*dJ1 - J1'*Tg1*l1
%     if ft+H*u_mu1(1)<0
%         Kxixi1 = w1*dJ1*Ge1'*C*Gc1;
%         Kt1 = Ke1 - Kdxi1*inv(Kxixi1)*Kxid1;
%     end
%     ri1 = B1'*sigma1*w1*dJ1+Kdxi1*(Kxixi1\renh1);
%     
%     ue2 = u(con2);
%     epsilon2 = B2*ue2 + Gc2*xi2;
%     sigma2 = C*epsilon2;
%     xi2
%     u_mu2 = J2*xi2
%     Tg2 = [max(0,ft+H*u_mu2(1));0];
%     renh2 = -Ge2'*sigma2*w2*dJ2 - J2'*Tg2*l2
%     if ft+H*u_mu2(1)<0
%         Kxixi2 = w2*dJ2*Ge2'*C*Gc2;
%         Kt2 = Ke2 - Kdxi2*inv(Kxixi2)*Kxid2;
%     end
%     ri2 = B2'*sigma2*w2*dJ2+Kdxi2*(Kxixi2\renh2);
%     
%     Kt = zeros(8);
%     Kt(con1,con1) = Kt(con1,con1) + Kt1;
%     Kt(con2,con2) = Kt(con2,con2) + Kt2;
%     
%     ri = zeros(8,1);
%     ri(con1) = ri(con1) + ri1;
%     ri(con2) = ri(con2) + ri2;
%     
%     Kt(free_dofs,free_dofs)
%     ri(free_dofs)
%     du = Kt(free_dofs,free_dofs)\(fext-ri(free_dofs));
%     u(free_dofs) = u(free_dofs) + du;
%     
%     Du(free_dofs) = du;
%     
%     due1 = Du(con1)
%     dxi1 = Kxixi1\(renh1-Kxid1*due1)
%     due2 = Du(con2);
%     dxi2 = Kxixi2\(renh2-Kxid2*due2)
%     
%     xi1 = xi1+dxi1;
%     xi2 = xi2+dxi2;
%     
%     conv_criteria = norm(du)/norm(u);
%     
% %     epsilon1
% %     epsilon2
%     i = i+1;
%     
%     cnt = cnt+1;
%     if(conv_criteria < tol) 
%         disp(['converged after ' ...
%                 num2str(cnt) ' loops']);
%         break;
%     end
% end
% u

end

